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Abstract 


For  a  data  assimilation  problem  for  a  limited-area  ocean  model,  the  boundary  data 
can  either  be  supplied  by  its  global  model  or  from  observations  interpolated  from 
conventional  observing  network.  In  this  study,  the  sensitivity  of  the  ocean  model  to 
variations  in  the  boundary  data  will  be  investigated.  Additionally,  we  will  be  looking 
at  how  predictions  of  wave  height  are  affected  by  variations  in  basin  depth.  Two 
models  are  used  to  carry  out  the  analyses;  a  large  scale  ocean  circulation  model  and 
a  combined  wave  refraction  diffraction  model.  The  large  scale  circulation  model  is 
used  to  implement  the  limited  area  ocean  model  and  results  show  that  the  accuracy 
and  availability  of  boundary  data  affects  the  model  predictions. 
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Chapter  1 


Introduction 

1.1  Background 

The  assimilation  of  random  data  into  ocean  models  has  not  been  extensively  looked  at. 
The  concern  of  most  researchers  is  that  the  models  axe  able  to  reproduce  the  behavior 
shown  by  collected  data.  A  data  assimilation  problem  involves  the  incorporation  of 
observed  data  into  the  problem,  usually  as  initial  conditions.  The  problem  is  then 
continually  updated  with  new  sets  of  observed  data.  Data  assimilation  also  involves 
the  preparation  of  observed  data  such  that  it  can  be  used  effectively  in  forecasting 
models.  It  is  an  effective  tool  to  have  when  available  data  is  limited.  For  a  data 
assimilation  problem  for  a  limited-area  ocean  model,  the  boundary  data  can  either 
be  supplied  by  its  global  model  or  from  observations  interpolated  from  conventional 
observing  network.  Either  way,  significant  uncertainty  is  generated  at  the  boundary 
of  the  model  representing  the  limited-area.  In  this  study,  the  sensitivity  of  the  ocean 
model  to  variations  in  the  boundary  data  will  be  investigated.  The  random  data  is 
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generated  from  a  coarse-grid  global  ocean  model  to  be  used  in  a  fine-grid  limited- 
area  model  with  open  boundaries.  Care  must  be  taken  in  how  the  boundary  data  is 
selected  as  it  has  been  shown  that  errors  in  the  initial  boundary  data  can  propagate 
into  the  interior  domain  significantly  affecting  the  stability  of  the  solution  [4]. 

Additionally,  the  sensitivity  of  wave  height  predictions  due  to  changes  in  bottom 
topography  are  also  investigated.  The  Karhunen-Loeve  expansion  is  used  to  generate 
random  data  to  be  used  for  the  bottom  topography  data.  The  bottom  is  thus  param¬ 
eterized,  as  a  stochastic  field  in  terms  of  a  few  random  amplitudes  associated  with 
depth.  Two  different  ocean  models  were  used  to  implement  the  objectives.  Ref/Dif 
1  developed  by  Kirby  and  Dalrymple  [6]  was  used  to  implement  the  bottom  topogra¬ 
phy.  The  general  ocean  circulation  model  developed  by  Gerald  Browning  at  NOAA1 
was  used  to  implement  the  open  boundary  conditions  for  a  limited  area  model.  The 
fundamental  difference  between  these  two  models  is  that  Ref/Dif  1  is  a  combined 
refraction  and  diffraction  model  that  models  the  propagation  of  water  waves  whereas 
the  NOAA  representation  models  the  general  ocean  circulation.  Hence,  in  Ref/Dif 
1  we  are  solving  for  the  wave  height  and  water  depth  and  in  the  general  circulation 
model,  for  the  current  velocity  and  pressure. 

1.2  Organization  of  this  Thesis 

This  thesis  is  organized  into  six  chapters.  Chapter  2  describes  how  a  general  ocean 
circulation  model  and  a  water  wave  model  are  developed.  Chapter  3  explains  the 
National  Oceanic  and  Atmospheric  Administration,  Boulder,  GO. 
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methods  used  for  the  sensitivity  analysis  in  both  models.  A  more  in-depth  look  into 
the  general  ocean  circulation  model  used  in  this  thesis  is  given  in  chapter  4.  Chapter 
4  also  describes  how  data  can  be  initialized  so  that  it  can  be  assimilated  into  the 
model.  The  results  of  the  sensitivity  analysis  on  the  limited  area  model  are  presented 
in  Chapter  4.  Chapter  5  describes  Ref/Dif  1.  It  also  includes  a  description  of  the 
analysis  carried  out  and  a  presentation  of  the  results.  Chapter  6  is  on  Conclusions  . 
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Chapter  2 


Model  Formulation 

2.1  Ocean  Circulation 

Ocean  currents  axe  driven  by  a  combination  of  the  sun,  wind  and  the  earth’s  rotation. 
The  most  obvious  interaction  that  drives  ocean  currents  is  that  between  the  wind  and 
the  sea  surface.  The  sun  also  drives  circulation  by  causing  variations  in  temperature 
and  salinity  which  result  in  a  density  change.  The  consequence  is  a  vertical  circulation 
pattern  called  the  thermohaline  circulation.  In  general,  warm  surface  currents  flow 
from  the  equator  to  the  poles  and  deep  water  cold  currents  flow  from  the  poles  to  the 
equator.  The  Coriolis  force,  which  is  a  result  of  the  earth  s  rotation,  is  zero  at  the 
equators  increasing  to  its  maximum  at  the  poles.  It  acts  at  right  angles  to  the  direction 
of  motion.  Hence  it  has  an  affect  on  the  direction  of  current  flow.  Furthermore, 
another  factor  affecting  the  direction  of  flow  is  land  boundaries.  The  combining 
effect  of  the  Coriolis  force  and  the  land  boundaries  causes  the  ocean  circulation  to 
be  in  a  circular  pattern.  The  circulation  model  developed  by  Munk  in  1950  closely 
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matches  the  actual  circulation  pattern  in  the  Atlantic  and  Pacific  Oceans.  The  basic 
equations  used  to  describe  the  circulation  are  briefly  reviewed  in  the  following  section. 

2.1.1  Governing  Equations 

In  the  derivation  of  the  equations  of  motion,  x  will  be  positive  in  the  eastward  direc¬ 
tion,  y  in  the  northward  and  z  is  positive  in  the  downward  direction.  The  velocity 
components  u,  v  and  w  are  in  the  x,  y,  and  z-direction  respectively.  The  equations 


Figure  2.1:  Coordinate  system  of  equations  of  motion 


used  to  describe  the  ocean  circulation  are  derived  from  the  conservation  of  momen¬ 
tum,  which  is  Newton’s  second  law  of  motion, 

F  =  ma.  (2.1) 


This  can  be  re-expressed  as 


accelaration  = 


force  per  unit  volume 
density 


(2,2) 
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Hence  the  equation  of  motion  in  the  x-direction  can  be  written  as 


du 

dt 


horizontal 

pressure 

gradient 

force  in 

x  —  direction 

Coriolis  force 
resulting  in 
motion  in 
x  —  direction 


other  forces  ] 
related  to 
motion  in 
x  —  direction  I 

/ 


(2.3) 


The  observation  that  the  Coriolis  force  acts  perpendicular  to  the  direction  of  flow  is 


represented  mathematically  as 


^  =  i(-^  +  p/v  +  Fe),  (2.4) 

dt  p  dx 

where  /  is  the  Coriolis  force, 1/  is  the  current  velocity  in  the  y-direction,  p  is  the 
pressure  and  p  is  the  density.  Fx  may  include  wind  stress,  friction  or  tidal  forcing. 
Similarly,  the  equation  of  motion  in  the  y-direction  is 

±  =  l(-±-pfu  +  F»).  (2.5) 

at  p  ay 

The  Coriolis  force  does  not  affect  flow  in  the  vertical  direction.  Hence  it  is  not  present 
in  the  vertical  equation  of  motion.  The  main  forces  in  the  vertical  direction  are  forces 
due  to  gravity  and  the  vertical  pressure  gradient.  Therefore  the  equation  of  motion 
in  the  z-direction  can  be  written  as 


^  =  I(_^_p£  +  F*).  (2.6) 

dt  p  dz 

As  this  problem  involves  an  advective  field,  where  a  change  in  the  field  is  a  result  of 
fluid  motion,  the  time  derivative  is  taken  to  be  the  material  derivative.  Hence, 


d  d  d  d  d 

Tt  =  m  +  uai  +  %  +  w8-z- 


(2.7) 


We  will  refer  to  the  above  equation  as  the  total  derivative  operator. 


6 


The  assumption  of  continuity  is  also  used.  The  continuity  equation  in  conjunction 
with  the  equations  of  motion  provide  additional  constraints  in  the  form  of  incompress¬ 


ibility. 


du  dv  9w_ 
dx  dy  dz 


(2.8) 


For  a  global  circulation  model  there  usually  is  land  on  the  lateral  boundaries. 
Hence  at  the  lateral  boundaries  the  velocity  normal  to  the  boundary  is  set  to  zero. 
Since  no  flow  is  allowed  through  the  sea  floor  at  the  sea  bottom,  the  vertical  velocity  is 
set  to  zero.  The  vertical  velocity  is  also  set  to  zero  at  the  sea  surface.  For  the  limited- 
area  model  the  same  boundary  conditions  will  apply  except  at  the  open  boundary 
where  the  flow  is  set  to  its  value  observed  in  the  global  model. 


2.2  Water  Waves 

The  wave  height,  water  depth  and  wave  length  or  wave  period  are  all  the  characteris¬ 
tics  needed  to  fully  describe  a  wave.  Figure  2.2  shows  a  two  dimensional  profile  of  a 
wave  with  these  characteristics  highlighted.  Other  wave  characteristics  such  as  water 
surface  profile,  the  forward  speed  of  the  wave  and  the  dynamic  pressure  field  in  the 
wave  can  be  theoretically  derived  from  the  aforementioned  characteristics. 

The  waves  that  are  of  interest  to  us  are  the  surface  gravity  waves.  With  these 
waves,  gravity  is  the  dominant  restoring  force.  Surface  tension  forces  are  neglected 
in  the  formulation  of  the  wave  equations.  This  is  because  surface  tension  forces  tend 
to  be  significant  only  when  the  waves  have  small  wavelengths,  for  example,  of  a 
couple  centimeters  and  since  surface  gravity  waves  have  wavelengths  on  the  order  of 
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Wave  Speed 


Length 


Depth 


Figure  2.2:  Two-dimensional  wave  profile 

a  couple  meters,  one  can  neglect  surface  tension  forces.  Some  of  the  processes  that 
waves  undergo  when  they  encounter  an  object  or  move  into  shallow  water  include 
shoaling,  refraction,  diffraction  and  dispersion.  Wave  refraction  is  the  change  in 
direction  of  propagation  which  occurs  when  isobaths,  which  are  lines  of  constant 
depth,  are  inclined  obliquely  to  wave  fronts  and  can  also  result  from  wave  current 
interactions.  Diffraction  is  the  change  of  wave  direction  due  to  changes  in  wave  height 
along  the  crest.  Dispersion  refers  to  the  separation  of  waves  due  to  their  different 
frequencies.  That  is,  waves  with  a  longer  wavelength  travel  faster  than  those  with  a 
shorter  wavelength.  These  three  phenomena  are  connected  through  the  dependence  of 
wave  energy  and  phase  velocity  on  depth  and  wave  height.  Wave  energy  is  dissipated 
as  a  result  of  boundary  shear  stresses  at  the  bottom  and  at  the  water-air  interface. 
A  porous  sand  bottom  will  allow  flow  to  occur  in  and  out  of  the  bottom.  This  can 
result  in  additional  energy  dissipation  if  the  water  is  shallow. 

In  the  development  of  wave  theories,  irrotational  flow  is  assumed.  This  is  a  valid 
assumption  since  viscous  forces  are  only  present  in  the  thin  boundary  layer  found 


8 


at  the  sea  surface  and  bottom,  despite  water  being  a  viscous  fluid.  The  governing 
equations  and  the  necessary  boundary  conditions  required  in  a  water  wave  model  are 
given  in  the  next  section. 


2.2.1  Governing  Equations 

Classical  linear  water  wave  theory  involves  the  following  assumptions  [2] 

1.  compressibility  and  viscosity  in  the  fluid  are  neglected  as  well  as  surface  tension 
at  the  free  surface; 

2.  the  fluid  motion  is  assumed  to  be  irrotational  and  can  be  described  by  a  velocity 
potential  <j>(x,  y,  z,  f); 

3.  nonlinear  terms  in  the  equations  of  motion  and  boundary  conditions  are  ne¬ 
glected. 


The  equations  resulting 


where 


from  the  above  assumptions  are, 

du  dv  dw 
dx  +  dy  dz 

(2.9) 

V  xU  =  0 

(2.10) 

U  =  V</> 

(2.11) 

then  V2(/>  =  0 

(2.12) 

U  =  (u(x,  y,  z,  t ),  v{x,  y,  z,  t ),  w(x,  y,  z,  t)). 

(2.13) 

Since  there  is  no  flow  normal  to  the  bottom,  the  bottom  boundary  condition  is 

u>  =  —  =  0  at  z  =  —d.  (2-14) 

dz 
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Two  boundary  conditions  can  exist  at  the  free  surface.  A  kinematic  boundary  condi¬ 
tion  and  a  dynamic  boundary  condition.  The  kinematic  boundary  condition  relates 
the  vertical  component  of  water  particle  velocity  at  the  surface  to  the  surface  position, 


dr)  dr) 
W  dt  +Udx 


at  z  =  r), 


(2.15) 


where  77  is  the  height  of  the  surface  above  the  still  water  level,  u  and  v  are  the  current 
velocity  in  the  x  and  y  -direction  respectively  and  t  is  time.  Nonlinear  terms  have 
been  neglected  according  to  linear  water  wave  theory.  Free  surfaces  cannot  support 
variations  in  pressure  across  the  air- water  interface  [8]  therefore  in  order  to  maintain  a 
uniform  pressure  along  the  interface  a  dynamic  boundary  condition  has  to  be  defined. 
The  dynamic  surface  boundary  condition  is  based  on  Bernoulli’s  equation  for  unsteady 
flow, 

^  +  !(t.2+t<2  +  i»2)  +  -+9Z  =  C(t),  (2.16) 

dt  2  p 

where  pv  is  the  pressure  at  the  free  surface,  p  is  the  water  density  and  </>  is  the 
velocity  potential.  The  current  velocity  components  are  represented  by  u,  v  and  w, 
g  the  acceleration  due  to  gravity  and  t  is  time.  The  pressure  at  the  free  surface  is 
constant  and  usually  taken  as  gauge  pressure,  pv  =  0.  C(t)  is  the  Bernoulli  term.  It 
is  constant  for  steady  flow.  Once  again  the  nonlinear  terms  have  been  left  out. 

When  the  depth  varies  with  x,  the  Laplace  equation  cannot  be  reduced  to  a  two 
dimensional  equation.  The  mild  slope  equation  was  developed  by  Berkhoff  (1972) 
to  handle  slowly  varying  bottom  topography.  The  mild  slope  equation  is  a  two  di¬ 
mensional  elliptical  partial  differential  equation  that  describes  the  complete  transfor¬ 
mation  of  small  amplitude  waves  including  refraction  and  diffraction.  The  Laplace 
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equation  is  vertically  integrated  to  give  the  following  equation, 


d  (nn  d<t>.  d  ,  d(f>  oCg 
Vx{CCTx)  +  Ty{CC^y)  +  a 


c 


=  0 


(2.17) 


where  C  is  the  wave  velocity,  Cg  is  the  wave  group  velocity  and  a  is  the  angular 
frequency.  The  wave  group  velocity  is  defined  as  the  overall  velocity  of  a  group  of 
ocean  waves  with  different  frequencies  and  wave  velocities  whereas,  the  wave  velocity 
is  the  velocity  of  a  single  wave.  The  wave  angular  frequency  ,  a  is  defined  as 


a2  =  gktaxihkd,  (2.18) 

where  k  is  the  wave  number,  d  is  the  water  depth  and  g  the  acceleration  due  to 
gravity. 
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Chapter  3 


Sensitivity  Analysis 


For  the  case  of  open  boundaries,  the  behavior  of  the  model  as  less  boundary  data  was 
available  was  investigated.  The  boundary  data  was  obtained  by  running  the  model 
over  a  closed  boundary  basin.  Simulations  were  then  run  with  varying  amounts  of 
boundary  data.  At  grid  points  where  no  boundary  data  was  available  a  linear  interpo¬ 
lation  was  carried  out.  Since  small  scale  disturbance  such  as  eddies  are  not  modeled 
by  the  ocean  circulation  model  the  linear  interpolation  technique  is  sufficient.  For  the 
case  of  modeling  the  bottom  topography  a  slightly  more  complex  methodology  was 
used.  It  involved  the  Karhunen-Loeve  expansion  and  Polynomial  Chaos  expansion  to 
generate  the  variability  in  the  bottom  topography  being  analyzed. 

3.1  Karhunen-Loeve  Expansion 

The  Karhunen-Loeve  expansion  of  a  stochastic  process  E(x ,  9)  is  based  on  the  spec¬ 
tral  expansion  of  its  covariance  function  Ree{x,v)  [10]  .where  x  and  y  are  spatial 
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coordinates  and  9  represents  the  random  nature  of  the  corresponding  quantity.  The 
expansion  is  defined  as, 


E(x, 6)  =  E(x)  +  £) 


(3.1) 


i=  1 


where  E(x )  is  the  mean  of  the  stochastic  process  and  &(0)  is  an  orthogonal  set  of 
random  variables.  The  eigenfunctions  and  eigenvalues  of  the  covariance  kernel  are 
represented  by  <&(x)  and  A,  respectfully.  The  eigenfunctions  are  mutually  orthogonal 
and  form  a  complete  set  spanning  the  function  space  to  which  E(x,9)  belongs  [10]. 
They  can  evaluated  as  the  solution  to  the  following  integral  equation 


[  REE{x-,y)(t>i{y)dy  =  \i<t>i{x),  (3.2) 

Jv 

where  D  is  the  domain  over  which  E(x,  9)  is  defined.  The  appeal  of  this  spectral 
representation  is  that  now  the  spatial  random  fluctuations  have  been  decomposed  into 
a  set  of  deterministic  functions  in  the  spatial  variables  multiplying  random  coefficients 
that  are  independent  of  these  variables.  When  a  process  is  closer  to  white  noise  more 
terms  are  needed  in  its  expansion  but  a  random  variable  can  be  represented  by  a 
single  term.  The  advantage  of  using  the  Karhunen-Loeve  expansion  is  that,  with 
only  a  few  terms  from  the  expansion  one  can  capture  most  of  the  uncertainty  in  the 
process,  since  it  can  be  assumed  that  material  properties  vary  smoothly. 


3.2  Polynomial  Chaos  Expansion 

Polynomial  Chaoses  are  polynomials  in  Gaussian  random  variables  that  describe  the 
functional  dependence  of  the  nodal  concentrations  c(9),  for  example  wave  height,  to 
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the  orthogonal  set  of  random  variables  &(0)  that  represents  the  material  stochasticity 
[10].  This  is  expressed  as 


OO  11 


c(6)  =  a0r0  +  £  +  £  £  ani2r2(4(0Ui2(0))  +  •  •  •  (3-3) 

i\  —  \  ii=l  i2=l 

where  rn(6x, . is  the  Polynomial  Chaos  of  order  n  in  the  Gaussian  variables 
(6i .  •  •  •  >  &„)•  The  pth  order  Polynomial  Chaos  can  be  expressed  as, 

rittm  t  ««_/ £L,(-ir&r(i, . .)re=ifi*(«)nr=r+i6,w  neven(34) 

rPte,w, 1  Eo U(-iriE.(j . iWnUi&(»)ns^.6,(«)  n odd  ^ 

where  n(-)  denotes  a  permutation  of  its  arguments  and  the  sum  is  over  all  such 
permutations  such  that  the  set  (4*n*  •  •  >Cv)  modified  by  the  permutation.  The 
Polynomial  Chaos  expansion,  when  truncated  after  the  pth  term  can  be  rewritten  as, 


c{6)  =  £c^(0) 

i= o 


(3.5) 


where  Cj  are  deterministic  coefficients  and  {^(0)}  is  a  basis  in  the  space  of  random 
variables.  This  basis  is  the  set  of  generalized  multidimensional  Hermite  polynomials. 
Due  to  the  orthogonality  of  the  Polynomial  Chaoses  the  deterministic  coefficient  can 
be  readily  obtained  by  multiplying  both  sides  of  equation  3.5  by  and  averaging. 
Hence,  the  deterministic  coefficients  are  determined  using  the  following  equation 


<  C^i  > 


Cj  <  *5  > 


(3.6) 


where  <  •  >  represents  the  mean.  The  denominator  of  the  equation  has  been  tabu¬ 
lated  [12]. 
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Chapter  4 

General  Circulation  Model 

4.1  Introduction 

This  section  refers  to  a  numerical  simulation  model  for  solving  the  general  circulation 
equations  with  specified  boundary  conditions.  The  model  describes  the  large  scale 
circulation  of  a  closed  basin.  At  the  present  it  is  run  over  a  square  basin  measuring 
2000  km  by  2000  km  by  2  km.  The  model  can  run  up  to  a  resolution  of  240  by  240 
by  24  grids  and  is  run  at  a  latitude  of  45  deg.  Hence  the  circulation  pattern  described 
by  the  model  closely  resembles  the  pattern  seen  in  the  North  Atlantic.  The  open 
boundaries  were  implemented  by  splitting  the  basin  in  half  such  that  one  has  outflow 
on  the  x  boundary.  The  shallow  water  equations  are  solved  in  the  model. 
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4.2  Shallow  Water  Equations 

The  model  uses  the  shallow  water  equations  to  describe  the  circulation.  The  shal¬ 
low  water  equations  are  a  system  of  hyperbolic  equations  describing  two  classes  of 
motion  with  different  time  scales.  The  “slow”  Rossby-type  (low  frequency)  motions 
and  “fast”  inertia-gravity  (high  frequency)  motions.  The  former  type  of  motions  are 
considered  to  be  significant,  as  they  contribute  the  most  to  the  energy  of  the  system. 
Several  methods  have  been  developed  to  filter  out  the  fast  waves.  In  initialization,  the 
underlying  equations  are  kept  but  the  initial  pressure  and  velocity  fields  are  adjusted 
so  as  to  control  the  high  frequency  motions  in  a  model.  One  initialization  method 
that  is  of  interest  is  the  bounded  derivative  method  developed  by  Kreiss  (1979).  The 
bounded  derivative  method  can  be  applied  to  both  global  and  limited-area  models. 
Other  initializations  schemes  such  as  the  nonlinear  normal  mode  initialization  pro¬ 
cedure  can  only  be  applied  to  a  limited  area  under  specific  boundary  conditions. 
Alternatively  the  equations  can  be  transformed  to  only  allow  low  frequency  motions. 

Browning  [4]  showed  that  only  the  large-scale  boundary  data  needs  to  be  correct 
in  order  to  prevent  the  boundary  error  from  propagating  into  the  interior  at  the 
speed  of  the  fast  waves.  The  process  of  specifying  data  along  the  boundary  of  the 
regional  model  can  be  viewed  as  interpolating  coarse  data  onto  a  fine  grid.  Given 
the  sensitivity  of  the  limited-area  model  to  certain  scales  of  boundary  information, 
it  is  expected  that  different  interpolation  schemes  may  result  in  different  modes  of 
error  propagation.  Ultimately,  it  may  be  desirable  especially  given  the  sparsity  of 
information  on  the  coarse  grid,  to  develop  good  stochastic  interpolation  schemes.  In 
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the  present  study,  however,  and  as  a  first  step  in  that  direction  the  significance  of  a 
variety  of  deterministic  interpolation  schemes  is  investigated. 

The  systems  of  equations  used  to  describe  shallow  water  behavior  are  usually 
derived  from  the  Eulerian  equations 


ds  _ 

jt-sw=F’ 

(4.1) 

du  _  i  ,  _ 

+  Po  Px-  fov  =  Fu 

(4.2) 

^  +  PolPy  +  fou  =  Fy 

(4.3) 

dw  i ,  .  „ 

-rr+Po  iPz  +  gs)  =  Fw 
at 

(4.4) 

+  c2p0(ux  +  vy  +  w2)  =  Fp 
at 

(4.5) 

d  d  d  d  d 

where  —  =  ^-  +  u—  +  v—  +  w— 
dt  dt  dx  ay  az 

(4.6) 

where  p0  is  the  mean  density,  p  the  pressure,  u,  v,  w  are  the  velocity  components,  f0 
ds  the  Coriolis  parameter.  F  is  the  forcing  function,  s  is  the  stratification  parameter, 
s  is  the  potential  density  and  c  is  the  speed  of  sound.  The  Coriolis  parameter  is  given 

by 

f  =  2Cl  sin  &o  +  —  cos  9q 
r 

where  Q  =  7.25  x  10_5s_1  is  the  earth’s  angular  speed,  r  is  the  earth’s  radius  and 
do  is  the  latitude  of  the  coordinate  origin.  For  our  case  9q  =  f  and  the  Coriolis 
parameter  is  assumed  to  be  constant.  The  potential  density  is  a  function  of  salinity 
and  potential  temperature.  The  potential  temperature  is  defined  as  the  temperature 
a  parcel  of  water  at  depth  would  attain  if  it  were  adiabatically  advected  to  the  surface. 
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The  potential  density  is  expressed  as 


s  =  999.84  +  6.79  x  1(T20  -  9.09  x  1(T302 

+  1.00  x  lO-403  -  1.12  x  1(T604  +  6.53  x  1(T905 
+  (0.824493  -  4.0899  x  1O_30  +  7.7438  x  lO-502  -  8.2467  x  1O_703 
+  5.3875  x  1O"904)5  +  (-5.72466  x  1(T3  +  1.0227  x  1O"40 
-  1.6546  x  1O_602)51-5  +  4.8314  x  lO-4-?2.  (4.7) 

In  the  above  equation  0  is  the  potential  temperature  in  deg  C  and  S  is  the  salinity  in 
practical  salinity  units.  By  making  different  assumptions,  variations  of  the  Eulerian 
system  can  be  obtained.  When  hydrostatic  and  incompressible  assumptions  are  made, 
the  primitive  system  of  equations  is  obtained.  When  the  gravity  and  sound  waves  are 
slowed  down  one  obtains  the  approximate  system. 

The  most  widely  used  systems  to  describe  the  shallow  water  equations  are 

1.  Eulerian  equations 

2.  Primitive  equations  -  derived  from  the  Eulerian  equations  by  applying  the  hy¬ 
drostatic  and  incompressible  assumptions.  A  drawback  of  this  system  is  that 
the  open  boundary  problem  is  always  ill  posed. 

3.  Reduced  system  -  derived  from  the  approximate  system  in  which  the  gravity 
and  sound  waves  have  been  slowed  down.  Hence  this  set  of  equations  only  allow 
the  low  frequency  motions  to  occur 

The  reduced  system  is  the  system  of  equations  used  in  the  NOAA  code.  The  scaled 
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reduced  system  is 


ds  ~  n 

— sw  =  0 

dt 

(4.8) 

~l(Px  -  fv)  =  Fu(y ,  z) 

(4.9) 

^  +  e_1(py  +  fu)  =  0 

(4.10) 

+  e~2(Pz  +  gs)  =  o 

(4.11) 

v*x  H"  Vy  +  ewz  0, 

(4.12) 

where  e  is  the  scaling  parameter  equal  to  the  Rossby  number  and  s  =  10“2. 

In  order  to  solve  the  shallow  water  equations,  they  first  need  to  be  scaled  due  to 
the  different  scales  of  motion  they  describe.  The  independent  variables  x,  y,  z  and 
time  are  renormalized  using  the  representative  length  and  time  scales.  The  dependent 
variables  are  scaled  by  their  representative  magnitudes.  The  shallow  water  equations 
are  solved  by  applying  the  total  time  derivative  to  the  continuity  equation.  An  elliptic 
equation  in  terms  of  the  pressure,  p  is  obtained  from  this  and  the  substitution  of  the 
equations  of  motion.  The  right  hand  side  of  the  elliptical  equation  is  a  function  of  s, 
u,  v,  w  and  their  first  derivatives  with  respect  to  x,  y  and  z. 

Pxx  +Pyy+Pzz  =  R(u,V,W,s),  (4.13) 

where 

agsz  df 

R(u,v,w,s)  =  f(vx  —  Uy)-\-2(uxVy+Ux‘wz-\-VyWz  —  uyvx  —  uzwx  —  vzWy)  -  ~u^' 

(4.14) 

The  forcing  function  does  not  appear  in  the  elliptical  equation  as  it  is  a  function  in  y 
and  z  only.  The  elliptical  equation  is  solved  using  an  elliptical  solver  from  FISHPAK. 
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The  velocity  components  and  potential  density  are  determined  numerically  using  a 
leapfrog  scheme, 

K+1  ~  _  K±±  ~  ffi-i  (4.15) 

2A  t  2A  h 

where  t  is  time,  k  is  the  property  being  solved  for  and  Ah  the  increment  of  the 
property. 

4.3  Boundary  Conditions 

The  global  model  is  an  initial  boundary  value  problem.  It  is  run  with  solid  boundaries 
on  all  four  sides.  At  the  boundaries  the  velocity  normal  to  the  boundary  is  set  to 
zero. 

u  =  0  @  x  =  west  boundary ,  x  =  east  boundary 
v  =  0  @  y  =  south  boundary,  y  =  north  boundary 

At  the  sea  surface  and  sea  bed  the  boundary  conditions  satisfy  the  condition  such 
that  no  water  passes  through  the  surface  and  sea  bed  respectively.  The  applied  stress 
is  zero  at  the  sea  bed.  Neumann  boundary  conditions  are  specified  for  the  global 
model 

=  fv  +  Fu  (4.16) 

dx 

±  =  -fu  (4.17) 

dy 

dv  =  Z^il  (4  18) 

dz  p 

For  the  limited-area  case,  Dirichlet  boundary  conditions  are  specified  on  the  open 
boundary.  These  conditions  are  obtained  from  post-processing  the  global  model  pre- 
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dictions. 


4.4  Bounded  Derivative  Method 

The  bounded  derivative  method  is  an  initialization  process  used  to  control  high  fre¬ 
quency  motions.  It  is  based  on  the  observation  that  the  solution,  which  varies  slowly 
with  respect  to  time,  of  a  hyperbolic  system  with  multiple  time  scales  must  have  a 
number  of  time  derivatives  on  the  order  of  the  slow  time  scale. 

The  bounded  derivative  method  can  be  applied  to  both  global  and  limited-area 
models.  The  principle  steps  used  to  determine  the  adjusted  initial  conditions  are  to, 

1.  define  the  characteristic  space  and  time  scales  of  the  motion  of  interest, 

2.  perform  a  scale  analysis  of  the  time  dependent  system  under  consideration  in 
order  to  identify  terms  that  can  contribute  to  large  time  derivatives, 

3.  constrain  those  terms  to  ensure  that  the  time  derivatives  are  of  the  low  time 
scale,  taking  into  account  the  boundary  conditions. 


4.5  Incorporating  Topography 


For  simplicity,  most  ocean,  models  are  developed  without  bottom  topography.  In 
order  to  incorporate  topography  into  the  approximate  systems  the  generalized  vertical 


coordinate  suggested  by  Kasahara  [11]  is  used.  It  has  the  advantage  of  keeping  the 
domain  of  integration  regular  and  of  allowing  considerable  freedom  in  the  placement 
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of  the  vertical  mesh  points.  The  independent  variables  are  transformed 


t  =  t' 

(4.19) 

X  =  X 

(4.20) 

y  =  y' 

(4.21) 

ii 

(4.22) 

The  velocity  components  are  obtained  by  applying  the  total  derivative  operator  to 

(4.20)-(4.22)  to  obtain 

u'  =  u 

(4,23) 

v'  —  V 

(4.24) 

w'  =  z'z(w  —  Zx>u'  —  Zy'v'). 

(4.25) 

In  the  transformed  coordinate  the  approximate  system  is 

ds  ~  n 

— - sw  =  0 

dt? 

(4.26) 

^7  +  Po  l(Px'  -  Zx'Z’zPz')  -fv  =  0 

(4.27) 

+  Po\Py>  -  Zy'Z'zpz>)  +  fu  =  0 

(4.28) 

^  +  a[polz'zpz'  +  giTPo^P  +  9s]  =  0 
at 

(4.29) 

^  +  lPo[ux'  -  zx>zzuz>  +  vx>  -  Zy>zzvz  +  z'zwz<)  -  gp0w  =  0, 
dt 

(4.30) 

where 

d  d  d  d  d  PP~lh  -  so 

M  =  di+UM+VW+W^'  S~  s° 

(4.31) 

The  system  is  valid  for  all  time  scales  because  the  unapproximated  forms  of  equations 
(4.27)  and  (4.28)  are  used.  In  the  large  scale  case  a  hydrostatic  approximation  is  used 
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to  replace  z'zpz>.  This  approximation  is  invalid  for  the  small  scale  case. 

To  develop  a  three  dimensional  flow,  the  untransformed  Eulerian  equations  are 
linearized  about  an  isothermal  ocean  at  rest.  This  results  in  the  linear  constant 
coefficient  system 

^--sw  =  0  (4.32) 

dt 

^  +  po  lpx  ~  fov  =  0  (4.33) 

^  +  Pq1Pv  +  /o«  =  0  (4.34) 

~  +  Po^Pz  +  gs  =  0  (4.35) 

at 

+  7Po(^x  +Vy+wz)  =  0  (4.36) 

at 


The  bounded  derivative  constraints  obtained  from  requiring  that  the  first  order 
time  derivatives  be  of  order  unity  is  that  the  flow  is  geostrophic,  hydrostatic  and 
totally  nondivergent.  Hence 


U  =  -'tpy 

(4.37) 

V  =  1px 

(4.38) 

(j)  =  fo-lp 

(4.39) 

-(SPo)~Vz> 

(4.40) 

where  'ip{x ,  y,  z)  is  a  given  streamfunction.  In  the  present  analysis,  the  bottom  topog 
raphy  in  the  circulation  model  is  assumed  constant  and  flat. 
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4.6  Sensitivity  Analysis 

A  global  model  with  a  coarse  grid  resolution  of  20  x  20  x  8  was  first  run  to  obtain 
the  necessary  boundary  data.  The  limited-area  model  was  a  basin  of  10  x  20  x  8, 
with  the  open  boundary  on  the  outflow  side.  Four  limited-area  cases  were  run.  In  the 
first  case,  all  the  boundary  data  from  the  global  model  was  passed  to  the  limited-area 
model.  In  the  second  case,  only  every  other  boundary  data  point  was  passed  to  the 
limited-area  model.  Linear  interpolation  was  used  to  approximate  the  data  at  the 
missed  points.  In  the  third  and  fourth  cases  only  five  and  three  boundary  data  points 
were  read  in  respectively.  Once  again  data  for  the  missed  grid  points  was  interpolated 
from  the  inputted  boundary  data. 

4.7  Results 

A  full  basin  model  of  grid  dimensions  20  x  20  x  8  was  run.  Using  the  contour  plots  of 
the  horizontal  velocity  u  as  an  indicator  the  implementation  of  the  open  boundaries 
seems  reasonable.  For  the  one  and  two  open  boundary  cases  the  behavior  agrees  at 
the  second  time  step.  Figures  4.1  -  4.3  show  contour  plots  of  a  full  basin  run,  a 
half-basin  run  with  the  open  boundary  on  the  eastern  boundary  and  a  three-quarters 
basin  run  with  open  boundaries  on  the  east  and  north  boundary.  From  Figures  4.2 
and  4.3  we  observe  an  exact  replication  of  the  global  circulation  in  a  limited-area.  A 
half-basin  is  used  to  analyses  the  behavior  of  the  open  boundary.  When  the  model 
is  left  to  run  for  longer  periods  we  see  an  exponential  increase  in  the  error.  The 
error  plot  for  our  four  open  boundary  conditions  is  shown  in  Figure  4.4.  The  solid 
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line  represents  data  taken  at  point  (6,16)  on  a  horizontal  plane.  The  dashed  line  is 
data  taken  at  point  (6,11).  As  mentioned  earlier,  only  the  boundary  data  needs  to 
be  accurate  in  order  to  prevent  error  from  destroying  the  solution.  However  we  see 
that  this  is  not  entirely  true.  In  the  error  plot  where  boundary  data  from  all  points 
were  passed  to  the  limited-area  model,  we  observe  error  entering  the  solution  at  an 
increasing  rate. 


Figure  4.1:  Contour  plot  of  the  horizontal  velocity  component,  u  over  the  entire  basin 
at  the  first  time  step 


BEST  AVAILABLE  COPY 
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2  day =  0  hr . =  0  min. -40  sec . =  0  k=  1 


canom  nw  «  M»ir~ 


Figure  4.2:  Contour  plot  of  the  horizontal  velocity  component,  u  over  half  the  basin 
with  one  open  boundary,  at  the  first  time  step 


2  day =  0  hr . =  0  min. =40  sec.=  0  K=  1 


CONTOUR  FROM  -4.5xl0~u  TO  4.5x10'**  BY  0.5x10“** 


Figure  4.3:  Contour  plot  of  the  horizontal  velocity  component,  u  over  three-quarters 
of  the  basin  with  two  open  boundaries  at  the  first  time  step 
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Chapter  5 

Water  Wave  Model 

5.1  Ref/Dif  1 

Ref/Dif  1  combines  both  the  effects  of  refraction  and  diffraction  explicitly,  making 
it  ideal  for  the  modeling  of  ocean  waves  in  regions  of  irregular  bottom  topography 
and  in  regions  where  diffraction  is  important  [6].  In  previous  models  in  order  to 
incorporate  diffraction,  refraction  was  suspended  in  the  regions  where  diffraction  was 
dominant.  This  method  was  inaccurate  but  it  allowed  for  the  inclusion  of  diffraction. 
Ray  tracing  techniques  were  used  for  the  analysis  of  refraction  of  water  waves.  These 
techniques  were  inaccurate  whenever  diffraction  effects  were  important,  since  they  do 
not  include  diffraction.  These  techniques  also  predicted  infinite  wave  heights  at  the 
point  of  intersecting  wave  rays  as  a  result  of  complex  bottom  topography.  Ref/Dif  1 
includes  the  following  effects  [6] 

1.  Parabolic  approximation: 
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(a)  Minimax  approximation 


2.  Wave  nonlinearity: 

(a)  Linear 

(b)  Composite  nonlinear 

(c)  Stokes  nonlinear 

3.  Wave  breaking 

4.  Absorbing  structures  and  shorelines: 

(a)  Thin  film  model  surrounded  by  a  natural  surf-zone 

5.  Energy  dissipation: 

(a)  Turbulent  bottom  friction  damping 

(b)  Porous  bottom  damping 

(c)  Laminar  boundary  layer  damping 

6.  Lateral  boundary  conditions: 

(a)  Reflective  condition 

(b)  Open  boundary  condition 

7.  Input  wave  field: 

(a)  Model  specification  of  monochromatic  or  directional  wave  field. 

(b)  Input  of  initial  row  of  data  from  disk  file. 


29 


8.  Output  wave  field: 


(a)  Standard  output 

(b)  Optional  storage  of  last  full  calculated  row  of  complex  amplitudes. 


5.2  Mild  slope  equation 


Very  few  solutions  exist  for  the  three  dimensional  problem  of  water  waves  propagat¬ 
ing  over  irregular  bathymetry.  This  problem  involves  complicated  nonlinear  bound¬ 
ary  conditions.  In  order  to  simplify  the  problem  it  was  reduced  to  a  problem  of  two 
horizontal  dimensions.  This  was  achieved  by  vertically  integrating  the  model.  Fur¬ 
thermore,  the  important  properties  of  linear  progressive  waves  could  be  predicted. 
The  linear  mild  slope  equation  is 


Tx(cc‘lb  +  +  a 


2%*=° 


where  C  is  the  wave  velocity,  Cg  is  the  group  velocity,  <fi  is  the  velocity  potential  and 
a  is  the  angular  frequency. 

The  model  equation  assumes  that  the  variations  in  the  bottom  occur  over  distances 
which  are  long  in  comparison  to  a  wave  length.  A  comparison  carried  out  by  Booij 
(1983)  found  that  the  mild  slope  model  was  accurate  for  bottom  slopes  up  to  1:3. 

The  model  equation  used  in  Ref/Dif  1  is  a  version  of  the  mild  slope  equation  that 
includes  the  influence  of  current.  The  equation  also  includes  the  minimax  approxi¬ 
mation,  which  is  the  minimization  of  the  maximum  error  for  a  rixed  number  of  terms 
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and  is  expressed  as, 


-  ib^y +  (7),  -  A2  {uU  (I). +  ru  (£)} 

+  iku>U(a0  —  1)  =  0,  (5.2) 


where 

fc  (fc(p  -  V *)), 

P  k 2  2k2 (p  -  U 2) 

Ai  =  Ox  —  &i 
A2  =  1  +  2ai  —  2 
A'  =  ai  -  61 1 

and  p  =  CCg,  U  is  the  mean  current  velocity  in  the  x  direction,  V  is  in  the  y 
direction,  k  is  the  average  wave  number  along  the  y  axis,  a;  is  a  dissipation  factor 
that  is  discussed  in  the  next  section  and  A  is  a  complex  amplitude  related  to  the 
water  surface  displacement  by 

77  =  Aei(fc*-at)  (5.3) 

In  the  above  equation  77  is  the  surface  displacement.  The  term  D  is  given  by, 

_  (cosh  4 kh  +  8-2  tanh2  kh)  ^ 

8  sinh4  kh 


31 


The  coefficients  ao,  cli  and  b\  depend  on  the  aperture  width  [6]  chosen  to  specify  the 
minimax  approximation.  However,  the  model  uses  the  Pade  approximation  based  on 


- - -  i  i 

the  following  coefficients 

ao  =  1 

(5.5) 

ai  =  —0.75 

(5.6) 

f>i  =  -0.25. 

(5.7) 

Since  the  model  follows  a  Stokes  expansion  it  is  only  valid  for  Stokes  waves.  Con¬ 
sequently  in  order  to  have  a  valid  model  in  shallow  water  a  heuristic  dispersion  re¬ 
lationship  developed  by  Hedges  (1976)  is  used.  The  final  model  has  a  dispersion 
relationship  that  is  a  smooth  patch  between  the  Hedges  form  and  Stokes  relationship. 


5.3  Energy  dissipation 


Ref/Dif  1  is  able  to  treat  energy  losses  due  to  bottom  friction,  surface  films  and  wave 
breaking.  The  type  of  energy  loss  determines  the  form  of  the  dissipation  factor,  u. 
For  turbulent  bottom  boundary  friction  the  dissipation  factor  is 

2akf\A\(l  -  i) 


u  = 


(5.8) 


37rsinh2A:/isinhA;/i 

The  Darcy- Weisbach  friction  factor,  /  is  set  to  0.01.  For  porous  bottom  damping  the 
dissipation  factor  is 

( i  „*\ 

(5.9) 


u> 


gkCp(l  —  ?) 


cosh2  kh 

where  Cp  is  the  coefficient  of  permeability  and  is  of  the  order  4.5  x  10~n  m2  and  g 
is  the  acceleration  due  to  gravity.  Laminar  boundary  layer  damping  occurs  both  at 
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the  water  surface  and  the  bottom.  It  is  a  result  of  viscosity.  Surface  boundary  layer 


damping  occurs  due  to  surface  films.  Surface  film  damping  is 

a\Ju]2a(l  -  i) 

U  tanh  kh 


(5.10) 


where  v  is  the  kinematic  viscosity.  The  term  under  the  square  root  sign  is  related  to 
the  thickness  of  the  boundary  layer  [6].  At  the  bottom  the  boundary  layer  damping 


is 


a 

u>  =  — 


The  wave  breaking  algorithm  is  always  active  in  the  model.  The  dissipation  due  to 
wave  breaking  is  given  by 


y/u/2a{l-i) 
sinh  kh 


(5.11) 


KCJ 1  -  (7 h/Hf) 


(5.12) 


where  K  and  7  are  arbitrary  constants  determined  to  equal  0.017  and  0.4  respectively. 


5.4  Bottom  Topography 

In  order  to  run  Ref/Dif  1  with  a  random  bottom,  the  physics  based  assumptions  used 
by  the  model  were  not  altered.  The  depth  data  in  the  output  data  file,  refdat.dat,  was 
changed  to  reflect  the  variable  bottom  topography.  The  Karhunen-Loeve  expansion 
was  used  to  represent  the  randomness  in  the  bottom  topography.  Ten  Gaussian 
variables,  &  were  used  to  represent  the  stochasticity  in  the  bottom  topography.  The 
elevation  of  the  bottom  was  consider  to  be  a  random  property  that  could  be  expanded 
as  follows, 

10 

k(x,e)  =  £fi(0)k,(x)  (5.13) 

1=0 
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where  k  is  the  elevation  of  the  bottom,  ki  represents  a  certain  scale  of  fluctuation  of 
bottom  elevation  and  represents  the  random  contribution  to  the  overall  property. 
The  10  eigenmodes  used  to  generate  the  bottom  topography  are  shown  in  Figures  7- 
17.  Various  solution  quantities  were  then  represented  in  terms  of  their  Polynomial 
Chaos  decomposition  as, 

cM,fO  =  X>OM)*i(0)  (5-14) 

i=0 

where  c<( x,  t)  are  deterministic  quantities  and  {'!/,(#)}  is  a  basis  in  the  space  of  random 
variables.  As  mentioned  earlier  in  our  discussion  of  the  Polynomial  Chaos  expansion, 
the  basis  is  taken  to  be  the  set  of  generalized  multidimensional  Hermite  Polynomials  in 
the  quantities  &(0).  Since  the  $  values  have  a  Gaussian  distribution  it  is  expected  that 
the  depth  of  the  model  will  also  follow  a  Gaussian  distribution.  The  Polynomial  chaos 
expansion  was  then  used  for  the  wave  height,  wave  angle  and  bottom  velocity  data. 
We  are  interested  in  seeing  the  type  of  distribution  these  quantities  will  now  have. 
Furthermore,  since  the  mild  slope  equation  is  accurate  for  slopes  not  greater  than  1:3, 
the  behavior  of  the  model  as  it  experiences  steeper  slopes  should  be  interpreted  with 
care. 

200  simulations  of  Ref/Dif  1  were  performed.  In  order  to  determine  the  coefficients 
Cj,  equation  5.14  was  divided  by  'k,  and  averaged 

<c^i>=Ci<^->  (5.15) 

The  right  hand  side  of  the  above  equation  can  be  expressed  as, 

<c#i>=  -t  £><>  (5.16) 

^  r=l 
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where  r  denotes  the  simulation  run.  After  employing  the  above  equation,  the  coeffi¬ 
cients  Ci  are  determined  using  equation  3.6. 
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5.5  Results 

Ref/Dif  1  was  run  to  simulate  the  behavior  of  waves  around  an  artificial  island.  200 
simulation  runs  were  performed.  The  island  is  circular  with  a  base  radius  of  400  ft. 
The  mean  water  depth  around  the  island  is  30  ft.  The  wave  conditions  modeled  were 
a  wave  height  of  28  ft  (see  figure  2)  and  a  wave  period  of  10  seconds.  On  input,  the  file 
containing  the  constant  depth  data  was  substituted  for  one  that  contained  variable 
depth  data.  Only  20  of  the  Hermite  Polynomials  were  used  to  generate  the  bottom 
topography.  Ref/Dif  1  produces  output  files  that  give  information  about  the  wave 
height,  the  water  depth  with  tide  correction,  wave  angle,  magnitude  of  the  bottom 
velocity  and  an  instantaneous  snapshot  of  the  water  surface.  Interest  lies  in  how  a 
randomly  generated  bottom  topography  will  affect  the  results  from  Ref/Dif  1.  Prior 
to  running  the  200  simulations,  the  model  of  the  island  was  run  with  a  flat  bottom 
topography.  The  results  from  this  run  are  shown  in  Figures  5.11  -  5.14.  In  Figures  5.12 
and  5.13  the  dashed  lines  are  contours  of  the  bottom  topography  and  the  solid  lines 
are  contours  of  wave  height  and  the  instantaneous  water  surface  respectively.  Due  to 
the  architecture  of  Ref/Dif  1,  we  assume  that  the  bottom  topography  to  the  left  of 
the  island  is  the  same  as  that  to  the  right  of  the  island.  Since  the  behavior  of  the 
water  waves  will  be  symmetric  only  half  the  island  is  shown  by  Ref/Dif  1.  The  waves 
are  incident  on  the  island  from  the  left  side  of  Figures  5.12-5.14. 

The  results  from  the  first  simulation  are  shown  in  Figures  5.15  -  5.17.  Figure  5.15 
shows  a  realization  of  the  bottom  topography.  A  plot  of  what  the  water  surface  looks 
like  is  given  by  Figure  5.16.  The  water  waves  are  incident  on  the  island  from  the  left 
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Figure  5.11:  Bottom  topography  for  flat  bottomed  model 

hand  side  of  Figure  5.18.  The  diffraction  of  the  water  fronts  as  they  encounter  the 
island  is  visible. 

The  probability  density  was  evaluated  for  the  wave  height,  water  depth,  wave 
angle  and  surface  (see  Figure  5.19).  The  probability  density  of  the  wave  angle  seems 
to  follow  a  Gaussian  distribution.  The  water  depth  shows  a  very  skewed  distribution 
as  expected  with  its  peak  occurring  around  31  ft.  The  density  plot  of  the  wave  height 
is  a  multimodal  distribution  with  three  peaks.  Its  highest  peak  occurs  at  12  ft. 

The  Polynomial  Chaos  expansion  was  truncated  after  the  21st  term  for  each  prop¬ 
erty.  The  first  term  in  the  expansion  represented  the  mean  of  the  property.  The  sub¬ 
sequent  summation  of  the  square  of  the  remaining  twenty  terms  yielded  the  variance 
of  the  property. 

Figures  5.24  -  5.27  are  plots  of  the  variance  for  the  properties. 
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Figure  5.18:  Rendering  of  water  surface  of  the  first  simulation 
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Figure  5.26:  Variance  of  wave  height 
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Figure  5.27:  Variance  of  surface 
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Chapter  6 
Conclusion 

For  the  limited-area  model,  error  was  observed  in  the  solution  even  when  all  the 
boundary  data  was  passed  to  the  limited-area  model.  This  hints  that  on  top  of  the 
need  for  accurate  boundary  data,  the  initial  data  might  need  to  be  initialized  by  the 
bounded  derivative  method,  for  example,  to  ensure  that  error  does  not  propagate  into 
the  solution. 

It  has  been  shown  that  variability  in  the  bottom  topography  does  not  affect  the 
results  obtained  from  Ref/Dif  1.  The  density  plots  of  the  wave  height  and  water  depth 
concur  with  the  mean  and  variance  plots  of  the  wave  height  and  bottom  topography. 
Furthermore,  a  comparison  of  the  mean  solution  to  the  deterministic  solution  shows 
that  the  random  bottom  topography  had  little  effect  on  the  Ref/Dif  1  results. 
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